%%capture

#pip install pretty-jupyter lxml_html_clean jinja-markdown setuptools scikit-learn scipy xgboost shap pandas numpy matplotlib seaborn probatus imblearn ipykernel
# jupyter nbconvert --execute --to html --template pj mynotebook.ipynb

%load_ext pretty_jupyter

import warnings
warnings.filterwarnings('ignore')

"""import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df_ct_rawhr = pd.read_pickle('clean_data/main_bigger3.pkl')
df_ct_rawhr['overallStatus'] = df_ct_rawhr['overallStatus'].replace({1:'TERMINATED', 0:'COMPLETED'})
df_ct_rawhr = df_ct_rawhr[df_ct_rawhr['is_Rare']==True]

# Function to calculate count and percentage
def calculate_counts_and_percentage(df, phase_name=None):
    if phase_name:
        phase_df = df[df[phase_name] == True]
    else:
        phase_df = df
    
    # Count of 'TERMINATED' (overallStatus = 0) and 'COMPLETED' (overallStatus = 1)
    counts = phase_df['overallStatus'].value_counts().reindex(['TERMINATED', 'COMPLETED'], fill_value=0)
    
    # Calculate percentage
    percentages = (counts / len(phase_df)) * 100
    
    return counts, percentages

# Function to calculate COVID trial counts
def calculate_covid_counts(df, phase_name=None):
    if phase_name:
        phase_df = df[df[phase_name] == True]
    else:
        phase_df = df
    
    covid_counts = phase_df[phase_df['is_covid_trial'] == 1]['overallStatus'].value_counts().reindex(['TERMINATED', 'COMPLETED'], fill_value=0)
    
    return covid_counts

# Filter data for 2019-2022
df_2019_2022 = df_ct_rawhr[(df_ct_rawhr['startDate'] >= 2019) & (df_ct_rawhr['startDate'] <= 2022)]

# Filter data for 2022 onwards
df_2022_onwards = df_ct_rawhr[df_ct_rawhr['startDate'] >= 2022]

# Calculate for all rows and phases for 2019-2022
all_counts_2019_2022, all_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022)
phase1_counts_2019_2022, phase1_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE1')
phase2_counts_2019_2022, phase2_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE2')
phase3_counts_2019_2022, phase3_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE3')

all_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022)
phase1_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE1')
phase2_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE2')
phase3_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE3')

# Calculate for all rows and phases for 2022 onwards
all_counts_2022_onwards, all_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards)
phase1_counts_2022_onwards, phase1_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE1')
phase2_counts_2022_onwards, phase2_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE2')
phase3_counts_2022_onwards, phase3_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE3')

all_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards)
phase1_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE1')
phase2_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE2')
phase3_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE3')

# Combine counts and percentages into DataFrames for easy plotting
plot_data_2019_2022 = pd.DataFrame({
    'Phase': ['All Rows', 'All Rows', 'PHASE1', 'PHASE1', 'PHASE2', 'PHASE2', 'PHASE3', 'PHASE3'],
    'Status': ['COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED'],
    'Count': [all_counts_2019_2022[1], all_counts_2019_2022[0], phase1_counts_2019_2022[1], phase1_counts_2019_2022[0],
              phase2_counts_2019_2022[1], phase2_counts_2019_2022[0], phase3_counts_2019_2022[1], phase3_counts_2019_2022[0]],
    'Percentage': [all_percentages_2019_2022[1], all_percentages_2019_2022[0], phase1_percentages_2019_2022[1], phase1_percentages_2019_2022[0],
                   phase2_percentages_2019_2022[1], phase2_percentages_2019_2022[0], phase3_percentages_2019_2022[1], phase3_percentages_2019_2022[0]],
    'Covid_Count': [all_covid_counts_2019_2022[1], all_covid_counts_2019_2022[0], phase1_covid_counts_2019_2022[1], phase1_covid_counts_2019_2022[0],
                    phase2_covid_counts_2019_2022[1], phase2_covid_counts_2019_2022[0], phase3_covid_counts_2019_2022[1], phase3_covid_counts_2019_2022[0]]
})

plot_data_2022_onwards = pd.DataFrame({
    'Phase': ['All Rows', 'All Rows', 'PHASE1', 'PHASE1', 'PHASE2', 'PHASE2', 'PHASE3', 'PHASE3'],
    'Status': ['COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED'],
    'Count': [all_counts_2022_onwards[1], all_counts_2022_onwards[0], phase1_counts_2022_onwards[1], phase1_counts_2022_onwards[0],
              phase2_counts_2022_onwards[1], phase2_counts_2022_onwards[0], phase3_counts_2022_onwards[1], phase3_counts_2022_onwards[0]],
    'Percentage': [all_percentages_2022_onwards[1], all_percentages_2022_onwards[0], phase1_percentages_2022_onwards[1], phase1_percentages_2022_onwards[0],
                   phase2_percentages_2022_onwards[1], phase2_percentages_2022_onwards[0], phase3_percentages_2022_onwards[1], phase3_percentages_2022_onwards[0]],
    'Covid_Count': [all_covid_counts_2022_onwards[1], all_covid_counts_2022_onwards[0], phase1_covid_counts_2022_onwards[1], phase1_covid_counts_2022_onwards[0],
                    phase2_covid_counts_2022_onwards[1], phase2_covid_counts_2022_onwards[0], phase3_covid_counts_2022_onwards[1], phase3_covid_counts_2022_onwards[0]]
})

# Plotting
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))

# Plot for 2019-2022
ax = axes[0]
x = np.arange(4)  # 4 phases (All Rows, PHASE1, PHASE2, PHASE3)
bar_width = 0.35  # Width of the bars

# Plot bars for COMPLETED
bars1 = ax.bar(x - bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Percentage'], 
               bar_width, label='COMPLETED')

# Plot bars for TERMINATED
bars2 = ax.bar(x + bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Percentage'], 
               bar_width, label='TERMINATED', color='orange')

# Plot bars for COVID trials inside the existing bars
covid_bars1 = ax.bar(x - bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Covid_Count'] / plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Count'] * plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Percentage'], 
                     bar_width, color='red', alpha=1, label='COVID Trials')

covid_bars2 = ax.bar(x + bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Covid_Count'] / plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Count'] * plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Percentage'], 
                     bar_width, color='red', alpha=1)

# Add text (counts) on top of bars
for bar, count, covid_count in zip(bars1, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Count'], plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

for bar, count, covid_count in zip(bars2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Count'], plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

# Set labels and title
ax.set_xlabel('Phases')
ax.set_ylabel('Percentage (%)')
ax.set_title('Percentage of COMPLETED and TERMINATED by Phase\nfor trials from 2019 to 2022')

ax.set_xticks(x)
ax.set_xticklabels(['All Rows', 'PHASE1', 'PHASE2', 'PHASE3'])
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.1), shadow=True, ncol=2)

# Plot for 2022 onwards
ax = axes[1]
x = np.arange(4)  # 4 phases (All Rows, PHASE1, PHASE2, PHASE3)
bar_width = 0.35  # Width of the bars

# Plot bars for COMPLETED
bars1 = ax.bar(x - bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Percentage'], 
               bar_width, label='COMPLETED')

# Plot bars for TERMINATED
bars2 = ax.bar(x + bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Percentage'], 
               bar_width, label='TERMINATED', color='orange')

# Plot bars for COVID trials inside the existing bars
covid_bars1 = ax.bar(x - bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Covid_Count'] / plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Count'] * plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Percentage'], 
                     bar_width, color='red', alpha=1, label='COVID Trials')

covid_bars2 = ax.bar(x + bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Covid_Count'] / plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Count'] * plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Percentage'], 
                     bar_width, color='red', alpha=1)

# Add text (counts) on top of bars
for bar, count, covid_count in zip(bars1, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Count'], plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

for bar, count, covid_count in zip(bars2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Count'], plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

# Set labels and title
ax.set_xlabel('Phases')
ax.set_ylabel('Percentage (%)')
ax.set_title('Percentage of COMPLETED and TERMINATED by Phase\nfor trials from 2022 onwards')

ax.set_xticks(x)
ax.set_xticklabels(['All Rows', 'PHASE1', 'PHASE2', 'PHASE3'])
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.1), shadow=True, ncol=2)

# Show plot
plt.tight_layout()
plt.show()"""

"""import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df_ct_rawhr = pd.read_pickle('clean_data/main_bigger3.pkl')
df_ct_rawhr['overallStatus'] = df_ct_rawhr['overallStatus'].replace({1:'TERMINATED', 0:'COMPLETED'})
df_ct_rawhr = df_ct_rawhr[(df_ct_rawhr['startDate'] >= 2019) & (df_ct_rawhr['startDate'] <= 2022)]

# Function to calculate count and percentage
def calculate_counts_and_percentage(df, phase_name=None):
    if phase_name:
        phase_df = df[df[phase_name] == True]
    else:
        phase_df = df
    
    # Count of 'TERMINATED' (overallStatus = 0) and 'COMPLETED' (overallStatus = 1)
    counts = phase_df['overallStatus'].value_counts().reindex(['TERMINATED', 'COMPLETED'], fill_value=0)
    
    # Calculate percentage
    percentages = (counts / len(phase_df)) * 100
    
    return counts, percentages

# Function to calculate COVID trial counts
def calculate_covid_counts(df, phase_name=None):
    if phase_name:
        phase_df = df[df[phase_name] == True]
    else:
        phase_df = df
    
    covid_counts = phase_df[phase_df['is_covid_trial'] == 1]['overallStatus'].value_counts().reindex(['TERMINATED', 'COMPLETED'], fill_value=0)
    
    return covid_counts

# Filter data for 2019-2022
df_2019_2022 = df_ct_rawhr[df_ct_rawhr['is_Rare']==True]

# Filter data for 2022 onwards
df_2022_onwards = df_ct_rawhr[df_ct_rawhr['is_Rare']==False]

# Calculate for all rows and phases for 2019-2022
all_counts_2019_2022, all_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022)
phase1_counts_2019_2022, phase1_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE1')
phase2_counts_2019_2022, phase2_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE2')
phase3_counts_2019_2022, phase3_percentages_2019_2022 = calculate_counts_and_percentage(df_2019_2022, 'PHASE3')

all_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022)
phase1_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE1')
phase2_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE2')
phase3_covid_counts_2019_2022 = calculate_covid_counts(df_2019_2022, 'PHASE3')

# Calculate for all rows and phases for 2022 onwards
all_counts_2022_onwards, all_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards)
phase1_counts_2022_onwards, phase1_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE1')
phase2_counts_2022_onwards, phase2_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE2')
phase3_counts_2022_onwards, phase3_percentages_2022_onwards = calculate_counts_and_percentage(df_2022_onwards, 'PHASE3')

all_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards)
phase1_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE1')
phase2_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE2')
phase3_covid_counts_2022_onwards = calculate_covid_counts(df_2022_onwards, 'PHASE3')

# Combine counts and percentages into DataFrames for easy plotting
plot_data_2019_2022 = pd.DataFrame({
    'Phase': ['All Rows', 'All Rows', 'PHASE1', 'PHASE1', 'PHASE2', 'PHASE2', 'PHASE3', 'PHASE3'],
    'Status': ['COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED'],
    'Count': [all_counts_2019_2022[1], all_counts_2019_2022[0], phase1_counts_2019_2022[1], phase1_counts_2019_2022[0],
              phase2_counts_2019_2022[1], phase2_counts_2019_2022[0], phase3_counts_2019_2022[1], phase3_counts_2019_2022[0]],
    'Percentage': [all_percentages_2019_2022[1], all_percentages_2019_2022[0], phase1_percentages_2019_2022[1], phase1_percentages_2019_2022[0],
                   phase2_percentages_2019_2022[1], phase2_percentages_2019_2022[0], phase3_percentages_2019_2022[1], phase3_percentages_2019_2022[0]],
    'Covid_Count': [all_covid_counts_2019_2022[1], all_covid_counts_2019_2022[0], phase1_covid_counts_2019_2022[1], phase1_covid_counts_2019_2022[0],
                    phase2_covid_counts_2019_2022[1], phase2_covid_counts_2019_2022[0], phase3_covid_counts_2019_2022[1], phase3_covid_counts_2019_2022[0]]
})

plot_data_2022_onwards = pd.DataFrame({
    'Phase': ['All Rows', 'All Rows', 'PHASE1', 'PHASE1', 'PHASE2', 'PHASE2', 'PHASE3', 'PHASE3'],
    'Status': ['COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED', 'COMPLETED', 'TERMINATED'],
    'Count': [all_counts_2022_onwards[1], all_counts_2022_onwards[0], phase1_counts_2022_onwards[1], phase1_counts_2022_onwards[0],
              phase2_counts_2022_onwards[1], phase2_counts_2022_onwards[0], phase3_counts_2022_onwards[1], phase3_counts_2022_onwards[0]],
    'Percentage': [all_percentages_2022_onwards[1], all_percentages_2022_onwards[0], phase1_percentages_2022_onwards[1], phase1_percentages_2022_onwards[0],
                   phase2_percentages_2022_onwards[1], phase2_percentages_2022_onwards[0], phase3_percentages_2022_onwards[1], phase3_percentages_2022_onwards[0]],
    'Covid_Count': [all_covid_counts_2022_onwards[1], all_covid_counts_2022_onwards[0], phase1_covid_counts_2022_onwards[1], phase1_covid_counts_2022_onwards[0],
                    phase2_covid_counts_2022_onwards[1], phase2_covid_counts_2022_onwards[0], phase3_covid_counts_2022_onwards[1], phase3_covid_counts_2022_onwards[0]]
})

# Plotting
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))

# Plot for 2019-2022
ax = axes[0]
x = np.arange(4)  # 4 phases (All Rows, PHASE1, PHASE2, PHASE3)
bar_width = 0.35  # Width of the bars

# Plot bars for COMPLETED
bars1 = ax.bar(x - bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Percentage'], 
               bar_width, label='COMPLETED')

# Plot bars for TERMINATED
bars2 = ax.bar(x + bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Percentage'], 
               bar_width, label='TERMINATED', color='orange')

# Plot bars for COVID trials inside the existing bars
covid_bars1 = ax.bar(x - bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Covid_Count'] / plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Count'] * plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Percentage'], 
                     bar_width, color='red', alpha=1, label='COVID Trials')

covid_bars2 = ax.bar(x + bar_width/2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Covid_Count'] / plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Count'] * plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Percentage'], 
                     bar_width, color='red', alpha=1)

# Add text (counts) on top of bars
for bar, count, covid_count in zip(bars1, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Count'], plot_data_2019_2022[plot_data_2019_2022['Status'] == 'COMPLETED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

for bar, count, covid_count in zip(bars2, plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Count'], plot_data_2019_2022[plot_data_2019_2022['Status'] == 'TERMINATED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

# Set labels and title
ax.set_xlabel('Phases')
ax.set_ylabel('Percentage (%)')
ax.set_title('Percentage of COMPLETED and TERMINATED by Phase\nfor Rare Diseases trials from 2019 to 2022')

ax.set_xticks(x)
ax.set_xticklabels(['All Rows', 'PHASE1', 'PHASE2', 'PHASE3'])
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.1), shadow=True, ncol=2)

# Plot for 2022 onwards
ax = axes[1]
x = np.arange(4)  # 4 phases (All Rows, PHASE1, PHASE2, PHASE3)
bar_width = 0.35  # Width of the bars

# Plot bars for COMPLETED
bars1 = ax.bar(x - bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Percentage'], 
               bar_width, label='COMPLETED')

# Plot bars for TERMINATED
bars2 = ax.bar(x + bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Percentage'], 
               bar_width, label='TERMINATED', color='orange')

# Plot bars for COVID trials inside the existing bars
covid_bars1 = ax.bar(x - bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Covid_Count'] / plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Count'] * plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Percentage'], 
                     bar_width, color='red', alpha=1, label='COVID Trials')

covid_bars2 = ax.bar(x + bar_width/2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Covid_Count'] / plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Count'] * plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Percentage'], 
                     bar_width, color='red', alpha=1)

# Add text (counts) on top of bars
for bar, count, covid_count in zip(bars1, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Count'], plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'COMPLETED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

for bar, count, covid_count in zip(bars2, plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Count'], plot_data_2022_onwards[plot_data_2022_onwards['Status'] == 'TERMINATED']['Covid_Count']):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 15, f"{int(count)}\n(COVID: {int(covid_count)})", ha='center', va='bottom')

# Set labels and title
ax.set_xlabel('Phases')
ax.set_ylabel('Percentage (%)')
ax.set_title('Percentage of COMPLETED and TERMINATED by Phase\nfor Non Rare Diseases trials from 2019 to 2022')

ax.set_xticks(x)
ax.set_xticklabels(['All Rows', 'PHASE1', 'PHASE2', 'PHASE3'])
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.1), shadow=True, ncol=2)

# Show plot
plt.tight_layout()
plt.show()"""

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Load the data
df_ct_rawhr = pd.read_pickle('clean_data/main_bigger3.pkl')
df_ct_rawhr['overallStatus'] = df_ct_rawhr['overallStatus'].replace({1:'TERMINATED', 0:'COMPLETED'})
df_ct_rawhr = df_ct_rawhr[(df_ct_rawhr['startDate'] >= 2019) & (df_ct_rawhr['startDate'] <= 2022)]
df_ct_rawhr = df_ct_rawhr[(df_ct_rawhr['PHASE2']==True)]

df_ct_rawhr_gen = pd.read_pickle('clean_data/main_bigger3.pkl')
df_ct_rawhr_gen['overallStatus'] = df_ct_rawhr_gen['overallStatus'].replace({1:'TERMINATED', 0:'COMPLETED'})
df_ct_rawhr_gen = df_ct_rawhr_gen[(df_ct_rawhr_gen['startDate'] >= 2022)]
df_ct_rawhr_gen = df_ct_rawhr_gen[(df_ct_rawhr_gen['PHASE2']==True)]

# Calculate value counts for overallStatus in both dataframes based on is_Rare
overall_status_rawhr = df_ct_rawhr.groupby(['is_Rare', 'overallStatus']).size().reset_index(name='count')
overall_status_rawhr_gen = df_ct_rawhr_gen.groupby(['is_Rare', 'overallStatus']).size().reset_index(name='count')

# Replace True/False with "Rare Diseases" and "Common Diseases"
overall_status_rawhr['is_Rare'] = overall_status_rawhr['is_Rare'].replace({True: 'Rare Diseases', False: 'Common Diseases'})
overall_status_rawhr_gen['is_Rare'] = overall_status_rawhr_gen['is_Rare'].replace({True: 'Rare Diseases', False: 'Common Diseases'})

# Ensure "Rare Diseases" comes first
overall_status_rawhr['is_Rare'] = pd.Categorical(overall_status_rawhr['is_Rare'], categories=['Rare Diseases', 'Common Diseases'])
overall_status_rawhr_gen['is_Rare'] = pd.Categorical(overall_status_rawhr_gen['is_Rare'], categories=['Rare Diseases', 'Common Diseases'])

# Plot the bar plots
plt.figure(figsize=(11, 5))

# Bar plot for df_ct_rawhr
plt.subplot(1, 2, 1)
ax1 = sns.barplot(data=overall_status_rawhr, x='is_Rare', y='count', hue='overallStatus', palette='viridis')
plt.title('Distribution of Phase 2 Rare Diseases and \nCommon Diseases Trials Status in the 2019 - 2022 Period') 
plt.xticks(rotation=0)
ax1.set(xlabel=None, ylabel=None)  # Hide x and y axis names

# Annotate bars with counts
for p in ax1.patches:
    if p.get_height() > 0:
        ax1.annotate(format(p.get_height(), '.0f'), 
                     (p.get_x() + p.get_width() / 2., p.get_height()), 
                     ha = 'center', va = 'center', 
                     xytext = (0, 5), 
                     textcoords = 'offset points')

# Bar plot for df_ct_rawhr_gen
plt.subplot(1, 2, 2)
ax2 = sns.barplot(data=overall_status_rawhr_gen, x='is_Rare', y='count', hue='overallStatus', palette='viridis')
plt.title('Distribution of Phase 2 Rare Diseases and \nCommon Diseases Trials Status in the Period beyond 2023') 
plt.xticks(rotation=0)
ax2.set(xlabel=None, ylabel=None)  # Hide x and y axis names

# Annotate bars with counts
for p in ax2.patches:
    if p.get_height() > 0:
        ax2.annotate(format(p.get_height(), '.0f'), 
                     (p.get_x() + p.get_width() / 2., p.get_height()), 
                     ha = 'center', va = 'center', 
                     xytext = (0, 5), 
                     textcoords = 'offset points')

plt.tight_layout()
plt.show()

import warnings
import pandas as pd
ct_raw = pd.read_pickle('clean_data/main_bigger3.pkl')

warnings.filterwarnings('ignore')

# print all the columns that are not boolean
non_bool_cols = ct_raw.columns[ct_raw.dtypes != 'bool']
# print the columns type


# Convert all other columns to boolean, except those in non_bool_cols
for col in ct_raw.columns:
    if col not in non_bool_cols:
        ct_raw[col] = ct_raw[col].astype(bool)
        
to_rmv=['enrolled_less_than_orig_enroll', 'less_than_75_orig_enroll','enrollmentCountSame','less_than_90_orig_enroll', 'less_than_50_orig_enroll', 'more_than_50_orig_enroll','more_than_75_orig_enroll','bb_mesh:Rare ']

ct_raw.drop(columns=to_rmv,inplace=True)

ct_raw = ct_raw[ct_raw['PHASE2']==True]
ct_rawr = ct_raw[ct_raw['is_Rare']==True]
ct_rawn = ct_raw[ct_raw['is_Rare']==False]

ct_rawr.drop(columns=['is_Rare','PHASE2'],inplace=True)
ct_rawn.drop(columns=['is_Rare','PHASE2'],inplace=True)









ct_fullr = ct_rawr.copy()
ct_fulln = ct_rawn.copy()



ct_rawr = ct_rawr[(ct_rawr['startDate'] >=2019)&(ct_rawr['startDate'] <=2022)]
ct_rawr.drop(columns=['startDate'],inplace=True)

ct_rawn = ct_rawn[(ct_rawn['startDate'] >=2019)&(ct_rawn['startDate'] <=2022)]
ct_rawn.drop(columns=['startDate'],inplace=True)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
import scipy.stats as stats
import numpy as np
from sklearn.preprocessing import FunctionTransformer
from sklearn.metrics import precision_recall_curve, roc_curve, auc, roc_auc_score
from collections import OrderedDict
from imblearn.over_sampling import RandomOverSampler, SMOTE
from sklearn.cluster import DBSCAN
from sklearn.compose import ColumnTransformer
from sklearn.datasets import fetch_openml
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectPercentile, chi2
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (accuracy_score, average_precision_score, classification_report, confusion_matrix, 
                             f1_score, precision_recall_curve, precision_score, recall_score, roc_auc_score, roc_curve)
from sklearn.model_selection import (GridSearchCV, RandomizedSearchCV, cross_val_score, learning_curve, 
                                     train_test_split)
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler,MinMaxScaler,RobustScaler,PowerTransformer,QuantileTransformer
from xgboost import XGBClassifier

warnings.filterwarnings('ignore')
pd.set_option('display.max_columns', None)

import pandas as pd
from sklearn.model_selection import train_test_split
from imblearn.over_sampling import RandomOverSampler
from sklearn.linear_model import LogisticRegression
from xgboost import XGBClassifier
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix, classification_report, precision_recall_curve, roc_curve, auc, roc_auc_score
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_validate
import seaborn as sns
import numpy as np

def LR_model(X, y):
    # Splitting the dataset into training and testing sets
    X_trainx, X_test, y_trainx, y_test = train_test_split(X, y, test_size=0.5, random_state=3)

    ros = RandomOverSampler(random_state=42)
    # Fit RandomOverSampler on the training data
    X_train, y_train = ros.fit_resample(X_trainx, y_trainx)
    feature_names = X_train.columns

    model = LogisticRegression(C=0.30001000000000005, penalty='l1', random_state=0, class_weight='balanced',solver='liblinear')
    model.fit(X_train, y_train)

    y_pred = model.predict(X_test)
    y_pred_prob = model.predict_proba(X_test)

    threshold = 0.53
    y_pred = (y_pred_prob[:, 1] >= threshold).astype(int)

    model_name = 'Logistic Regression'
    accuracy = accuracy_score(y_test, y_pred)
    xprecision = precision_score(y_test, y_pred)
    xrecall = recall_score(y_test, y_pred)
    f1 = f1_score(y_test, y_pred)
    conf_matrix = confusion_matrix(y_test, y_pred)

    print(f"Accuracy: {accuracy:.3f}")
    print(f"Precision: {xprecision:.3f}")
    print(f"Recall: {xrecall:.3f}")
    print(f"F1 Score: {f1:.3f}")
    print("Confusion Matrix:")

    plt.figure(figsize=(18, 6))
    plt.subplot(1, 3, 1)
    sns.heatmap(conf_matrix, annot=True, fmt='g', cmap='Blues', xticklabels=['Predicted Negative', 'Predicted Positive'], yticklabels=['Actual Negative', 'Actual Positive'])
    plt.xlabel('Predicted labels')
    plt.ylabel('True labels')
    plt.title('Confusion Matrix for ' + model_name)

    print("Classification Report:")
    print(classification_report(y_test, y_pred))

    precision, recall, _ = precision_recall_curve(y_test, y_pred_prob[:, 1])
    pr_auc = auc(recall, precision)

    fpr, tpr, _ = roc_curve(y_test, y_pred_prob[:, 1])
    roc_auc = auc(fpr, tpr)

    print(f"PR-AUC: {pr_auc:.3f}")
    print(f"ROC-AUC: {roc_auc:.3f}")

    plt.subplot(1, 3, 2)
    plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (area = {roc_auc:.2f})')
    plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver Operating Characteristic Curve for ' + model_name)
    plt.legend(loc="lower right")

    plt.subplot(1, 3, 3)
    plt.plot(recall, precision, color='blue', lw=2, label=f'PR curve (area = {pr_auc:.2f})')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.title('Precision-Recall Curve for ' + model_name)
    plt.legend(loc="lower left")

    plt.show()

  
    # Define the scoring metrics
    scoring = {
        'accuracy': 'accuracy',
        'log_loss': 'neg_log_loss',
        'auc': 'roc_auc',
        'f1': 'f1',
        'pr_auc': 'average_precision',
        'recall': 'recall',
        'precision' : 'precision'
    }

    # Perform cross-validation
    results = cross_validate(model, X, y, cv=10, scoring=list(scoring.values()), return_train_score=False)

    print('K-fold cross-validation results:')
    metrics_list = []

    # Iterate over the scoring metrics
    for sc in range(len(scoring)):
        metric_name = list(scoring.keys())[sc]
        mean_score = -results['test_%s' % list(scoring.values())[sc]].mean() if list(scoring.values())[sc] == 'neg_log_loss' else results['test_%s' % list(scoring.values())[sc]].mean()
        std_score = results['test_%s' % list(scoring.values())[sc]].std()
        metrics_list.append({'Metric': metric_name, 'Mean Score': mean_score, 'Standard Deviation': std_score})
        print(model.__class__.__name__ + " average %s: %.3f ± %.2f" % (metric_name, mean_score, std_score))

    # Convert the list of metrics to a DataFrame
    cv_df = pd.DataFrame(metrics_list)

    # get a prediction df
    pred_df = pd.DataFrame(y_test)
    pred_df.reset_index(drop=False, inplace=True)

    X_train = X_train.astype(int)
    X_test = X_test.astype(int)
    # append values of y_pred to y_test
    pred_df['predicted'] = y_pred
    results = {
        'pred_df': pred_df,
        'X_train': X_train,
        'X_test': X_test,
        'y_train': y_train,
        'y_test': y_test,
        'accuracy': accuracy,
        'precision': xprecision,
        'recall': xrecall,
        'f1': f1,
        'pr_auc': pr_auc,
        'roc_auc': roc_auc,
        'model': model,
        'CV_scoring': cv_df
    }
    return results

def XGB_model(X,y):
    # Splitting the dataset into training and testing sets
    X_trainx, X_test, y_trainx, y_test = train_test_split(X, y, test_size=0.5, random_state=3)

    ros = RandomOverSampler(random_state=42)
    # Fit RandomOverSampler on the training data
    X_train, y_train = ros.fit_resample(X_trainx, y_trainx)
    feature_names = X_train.columns


    model = XGBClassifier(random_state=42, use_label_encoder=True, eval_metric='logloss', class_weight='balanced',
                          learning_rate= 0.1, max_depth= 5, n_estimators= 300, subsample= 0.8)

    model.fit(X_train, y_train)

    y_pred = model.predict(X_test)
    y_pred_prob = model.predict_proba(X_test)


    threshold = 0.6
    y_pred= (y_pred_prob[:, 1] >= threshold).astype(int)

    model_name = 'XGBoost'
    accuracy = accuracy_score(y_test, y_pred)
    xprecision = precision_score(y_test, y_pred)
    xrecall = recall_score(y_test, y_pred)
    f1 = f1_score(y_test, y_pred)
    conf_matrix = confusion_matrix(y_test, y_pred)

    print(f"Accuracy: {accuracy:.3f}")
    print(f"Precision: {xprecision:.3f}")

    print(f"Recall: {xrecall:.3f}")
    print(f"F1 Score: {f1:.3f}")
    print("Confusion Matrix:")

    # Plotting the confusion matrix
    plt.figure(figsize=(18,6))
    plt.subplot(1, 3, 1)  # 1 row, 2 columns, 1st subplot
    sns.heatmap(conf_matrix, annot=True, fmt='g', cmap='Blues', xticklabels=['Predicted Negative', 'Predicted Positive'], yticklabels=['Actual Negative', 'Actual Positive'])
    plt.xlabel('Predicted labels')
    plt.ylabel('True labels')
    plt.title('Confusion Matrix for '+ model_name)

    # Detailed classification report
    print("Classification Report:")
    print(classification_report(y_test, y_pred))

    # Calculate precision-recall pairs for different probability thresholds
    precision, recall, _ = precision_recall_curve(y_test, y_pred_prob[:, 1])
    pr_auc = auc(recall, precision)

    # Calculate the ROC curve
    fpr, tpr, _ = roc_curve(y_test, y_pred_prob[:, 1])
    roc_auc = auc(fpr, tpr)


    print(f"PR-AUC: {pr_auc:.3f}")
    print(f"ROC-AUC: {roc_auc:.3f}")



    # Plotting ROC Curve in the first subplot
    plt.subplot(1, 3, 2)  # 1 row, 2 columns, 1st subplot
    plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (area = {roc_auc:.2f})')
    plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver Operating Characteristic Curve for '+ model_name)
    plt.legend(loc="lower right")

    # Plotting Precision-Recall Curve in the second subplot
    plt.subplot(1, 3, 3)  # 1 row, 2 columns, 2nd subplot
    plt.plot(recall, precision, color='blue', lw=2, label=f'PR curve (area = {pr_auc:.2f})')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.title('Precision-Recall Curve for '+ model_name)
    plt.legend(loc="lower left")

    plt.show()



    # Define the scoring metrics
    scoring = {
        'accuracy': 'accuracy',
        'log_loss': 'neg_log_loss',
        'auc': 'roc_auc',
        'f1': 'f1',
        'pr_auc': 'average_precision',
        'recall': 'recall',
        'precision' : 'precision'
    }

    # Perform cross-validation
    results = cross_validate(model, X, y, cv=10, scoring=list(scoring.values()), return_train_score=False)

    print('K-fold cross-validation results:')
    metrics_list = []

    # Iterate over the scoring metrics
    for sc in range(len(scoring)):
        metric_name = list(scoring.keys())[sc]
        mean_score = -results['test_%s' % list(scoring.values())[sc]].mean() if list(scoring.values())[sc] == 'neg_log_loss' else results['test_%s' % list(scoring.values())[sc]].mean()
        std_score = results['test_%s' % list(scoring.values())[sc]].std()
        metrics_list.append({'Metric': metric_name, 'Mean Score': mean_score, 'Standard Deviation': std_score})
        print(model.__class__.__name__ + " average %s: %.3f ± %.2f" % (metric_name, mean_score, std_score))

    # Convert the list of metrics to a DataFrame
    cv_df = pd.DataFrame(metrics_list)

    # get a prediction df
    pred_df = pd.DataFrame(y_test)
    pred_df.reset_index(drop=False, inplace=True)

    X_train = X_train.astype(int)
    X_test = X_test.astype(int)
    # append values of y_pred to y_test
    pred_df['predicted'] = y_pred
    results = {
        'pred_df': pred_df,
        'X_train': X_train,
        'X_test': X_test,
        'y_train': y_train,
        'y_test': y_test,
        'accuracy': accuracy,
        'precision': xprecision,
        'recall': xrecall,
        'f1': f1,
        'pr_auc': pr_auc,
        'roc_auc': roc_auc,
        'model': model,
        'CV_scoring': cv_df
    }
    return results

import numpy as np
import pandas as pd
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, precision_recall_curve, auc, confusion_matrix
from sklearn.model_selection import cross_validate

def plot_validation(dict, ct_full):
    # Filter dataframe for startDate 2023 and 2024
    df_filtered = ct_full[ct_full['startDate'].isin([2023, 2024])]
    df_filtered.drop('startDate', axis=1, inplace=True)
    
    X = df_filtered.drop('overallStatus', axis=1)
    y = df_filtered['overallStatus']
    
    model = dict['model']
    metrics = ['precision','f1','auc','pr_auc','accuracy', 'recall', 'log_loss']
    # Define the scoring metrics
    scoring = {
        'precision' : 'precision',
        'f1': 'f1',
        'auc': 'roc_auc',
        'pr_auc': 'average_precision',
        'accuracy': 'accuracy',
        'recall': 'recall',
        'log_loss': 'neg_log_loss'
    }

    # Perform cross-validation
    results = cross_validate(model, X, y, cv=10, scoring=list(scoring.values()), return_train_score=False)

    print('K-fold cross-validation results:')
    metrics_list = []

    # Iterate over the scoring metrics
    for sc in range(len(scoring)):
        metric_name = list(scoring.keys())[sc]
        mean_score = -results['test_%s' % list(scoring.values())[sc]].mean() if list(scoring.values())[sc] == 'neg_log_loss' else results['test_%s' % list(scoring.values())[sc]].mean()
        std_score = results['test_%s' % list(scoring.values())[sc]].std()
        metrics_list.append({'Metric': metric_name, 'Mean Score': mean_score, 'Standard Deviation': std_score})
        print(model.__class__.__name__ + " average %s: %.3f ± %.2f" % (metric_name, mean_score, std_score))
    
    # Create and print the DataFrame with metrics
    metrics_df = pd.DataFrame(metrics_list)
    print(metrics_df)

## print LR coefficients
def get_coeff(model,feature_names,col_key_name):

    coefficients = model.coef_[0]  # Coefficients for the features

    coef_df = pd.DataFrame({'Feature': feature_names, 'Coefficient': coefficients})
    coef_df['Absolute Coefficient'] = coef_df['Coefficient'].abs()
    coef_df = coef_df.sort_values(by='Absolute Coefficient', ascending=False)
    coef_df['Absolute_Coefficient'] = coef_df['Coefficient'].abs()
    
    if col_key_name == 'top30':
        # Sort the DataFrame based on the absolute values in descending order and take the top 30
        top_30 = coef_df.sort_values(by='Absolute_Coefficient', ascending=False).head(30)

        # Plotting with feature names as y-axis labels
        plt.figure(figsize=(10, 8))
        top_30.plot(kind='barh', x='Feature', y='Coefficient', color='skyblue')
        plt.xlabel('Coefficient Value')
        plt.title('Top 30 Absolute Coefficient Values')
        plt.gca().invert_yaxis()  # To display the highest value at the top
        plt.show()

    
    # Existing code for creating lists for each prefix
    prefixes = ["continent", "country", "condition_keyword:", "mesh","bb_mesh", "org_name", "lead_sponsor","float"]
    columns_by_prefix = {prefix: [] for prefix in prefixes}

    # Populate the lists with column names that start with the respective prefixes
    for col in ct_raw.columns:
        matched = False
        for prefix in prefixes:
            if col.startswith(prefix):
                columns_by_prefix[prefix].append(col)
                matched = True
                break
        if not matched and col != 'overallStatus' and col != 'startDate':
            # If the column doesn't match any prefix, add it to a 'remaining' list
            if 'remaining' not in columns_by_prefix:
                columns_by_prefix['remaining'] = []
            columns_by_prefix['remaining'].append(col)

    # Create a new ordered dictionary with 'remaining' first, if it exists
    ordered_columns_by_prefix = OrderedDict()
    if 'remaining' in columns_by_prefix:
        ordered_columns_by_prefix['remaining'] = columns_by_prefix['remaining']
    for prefix in prefixes:
        ordered_columns_by_prefix[prefix] = columns_by_prefix[prefix]

    # Use ordered_columns_by_prefix for ordered access

    # Step 1: Define a mapping of old keys to new keys
    key_mapping = {
        'remaining': 'main',
        'continent': 'Continents compared to Western Europe',
        'country': 'Countries compared to France',
    }

    # Step 2: Iterate over the old dictionary to create a new one with the new keys
    new_ordered_columns_by_prefix = {key_mapping.get(old_key, old_key): value for old_key, value in ordered_columns_by_prefix.items()}

    # Now, new_ordered_columns_by_prefix contains the updated keys
    for name, prefixes in new_ordered_columns_by_prefix.items():
        if name == col_key_name:
            # Filter coef_df based on whether 'Feature' starts with any of the specified prefixes
            relevant_columns = [col for col in coef_df['Feature'] if any(col.startswith(prefix) for prefix in prefixes)]
            
            # Now, filter coef_df to only include rows where 'Feature' is in relevant_columns
            filtered_coef_df = coef_df[coef_df['Feature'].isin(relevant_columns)]
            filtered_coef_df = filtered_coef_df[filtered_coef_df['Absolute Coefficient'] > 0]
            filtered_coef_df = filtered_coef_df.sort_values(by='Coefficient')

            filtered_coef_df_positive = filtered_coef_df[filtered_coef_df['Coefficient'] > 0]
            filtered_coef_df_negative = filtered_coef_df[filtered_coef_df['Coefficient'] < 0]

            # Check if filtered_coef_df_positive is not empty
            if not filtered_coef_df_positive.empty:
                plt.figure(figsize=(10, 12))
                plt.rcParams.update({'font.size': 16})
                plt.barh(filtered_coef_df_positive['Feature'], filtered_coef_df_positive['Coefficient'])
                plt.xlabel('Coefficient Value')
                plt.ylabel('Feature')
                plt.title(f'{name} Feature Coefficients that increase probability of early Termination')
                plt.grid(True)
                plt.show()

            # Check if filtered_coef_df_negative is not empty
            if not filtered_coef_df_negative.empty:
                plt.figure(figsize=(10, 12))
                plt.rcParams.update({'font.size': 16})
                plt.barh(filtered_coef_df_negative['Feature'], filtered_coef_df_negative['Coefficient'])
                plt.xlabel('Coefficient Value')
                plt.ylabel('Feature')
                plt.title(f'{name} Feature Coefficients that decrease probability of early Termination')
                plt.grid(True)
                plt.show()
    if col_key_name == "float":
            # Filter coef_df based on whether 'Feature' starts with any of the specified prefixes
            
            # Now, filter coef_df to only include rows where 'Feature' is in relevant_columns
            filtered_coef_df = coef_df[coef_df['Feature'].isin(non_bool_cols)]
            filtered_coef_df = filtered_coef_df[filtered_coef_df['Absolute Coefficient'] > 0]
            filtered_coef_df = filtered_coef_df.sort_values(by='Coefficient')

            # Check if filtered_coef_df_positive is not empty
            if not filtered_coef_df.empty:
                plt.figure(figsize=(10, 12))
                plt.rcParams.update({'font.size': 16})
                plt.barh(filtered_coef_df['Feature'], filtered_coef_df['Coefficient'])
                plt.xlabel('Coefficient Value')
                plt.ylabel('Feature')
                plt.title(f'{name} Feature Coefficients that increase probability of early Termination')
                plt.grid(True)
                plt.show()
    if col_key_name == "masking":
        # Filter coef_df based on whether 'Feature' starts with any of the specified prefixes
        
        # Now, filter coef_df to only include rows where 'Feature' is in relevant_columns
        filtered_coef_df = coef_df[coef_df['Feature'].isin(masking_col)]
        filtered_coef_df = filtered_coef_df[filtered_coef_df['Absolute Coefficient'] > 0]
        filtered_coef_df = filtered_coef_df.sort_values(by='Coefficient')

        # Check if filtered_coef_df_positive is not empty
        if not filtered_coef_df.empty:
            plt.figure(figsize=(10, 12))
            plt.rcParams.update({'font.size': 16})
            plt.barh(filtered_coef_df['Feature'], filtered_coef_df['Coefficient'])
            plt.xlabel('Coefficient Value')
            plt.ylabel('Feature')
            plt.title(f'{name} Feature Coefficients that increase probability of early Termination')
            plt.grid(True)
            plt.show()
          

# list of all columns startin with the prefixe "continent"
continent_col = [col for col in ct_raw.columns if col.startswith('continent')]
country_col = [col for col in ct_raw.columns if col.startswith('country')]
mesh_col = [col for col in ct_raw.columns if col.startswith('mesh')]
org_col = [col for col in ct_raw.columns if col.startswith('org_name')]
lead_col = [col for col in ct_raw.columns if col.startswith('lead_sponsor_')]
bb_mesh = [col for col in ct_raw.columns if col.startswith('bb_mesh')]
### ordinal columns

phases_col = ['PHASE1','PHASE2','PHASE3']
enrollment_col_r=['original_enrollment_count','enrolled_less_than_orig_enroll', 'less_than_90_orig_enroll', 'less_than_75_orig_enroll', 'less_than_50_orig_enroll', 'more_than_50_orig_enroll']
exp_drug_col =['experimental_only','exp_against_placebo','exp_against_active_comp']
collaborators_col = [col for col in ct_raw.columns if col.startswith('collaborators =_')]
centers_col = [col for col in ct_raw.columns if col.startswith('total_centers_count =_')]
enrollment_col = [col for col in ct_raw.columns if col.startswith('enrollment_count_cat_')]
entity_col = [col for col in ct_raw.columns if any(col.startswith(prefix) for prefix in ['sponsor_class_', 'org_study_class'])]

intervention_col = [col for col in ct_raw.columns if col.startswith('interventionModel_')]

primarypurpose_col = [col for col in ct_raw.columns if any(col.startswith(prefix) for prefix in ['primaryPurpose', 'intervention_type_'])]
masking_1 = [col for col in ct_raw.columns if any(col.startswith(prefix) for prefix in ['mask_type_'])]
masking_2 = [col for col in ct_raw.columns if any(col.endswith(suffix) for suffix in ['_masked'])]
masking_col = masking_1+masking_2

excluded_cols = masking_col+phases_col +bb_mesh+ collaborators_col + centers_col + enrollment_col + continent_col + country_col + mesh_col + org_col + lead_col + entity_col + intervention_col + primarypurpose_col
main_col = [col for col in ct_raw.columns[1:] if col not in excluded_cols]
categ_col = collaborators_col+centers_col+enrollment_col

# drop columns from categ_col
#ct_raw.drop(columns=categ_col, inplace=True)

# List of DataFrames
dfs = [ct_rawr, ct_fullr, ct_rawn, ct_fulln]

# Function to drop columns with all the same values, except for 'overallStatus'
def drop_constant_columns(df):
    # Identify columns with more than one unique value or the column 'overallStatus'
    cols_to_keep = df.nunique() > 1
    cols_to_keep['overallStatus'] = True
    return df.loc[:, cols_to_keep]

# Apply the function to each DataFrame
ct_rawr, ct_fullr, ct_rawn, ct_fulln = [drop_constant_columns(df) for df in dfs]

# Assuming 'overallStatus' is your target variable and the rest are features

X = ct_rawr.drop('overallStatus', axis=1)
y = ct_rawr['overallStatus']
logreg_noise = LR_model(X,y)

# Assuming 'overallStatus' is your target variable and the rest are features
X1 = ct_rawn.drop('overallStatus', axis=1)
y1 = ct_rawn['overallStatus']
logreg_noise1 = LR_model(X1,y1)

# https://ing-bank.github.io/probatus/api/model_interpret.html
# pip install probatus
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from probatus.interpret import ShapModelInterpreter
import numpy as np
import pandas as pd
import shap
import matplotlib.pyplot as plt

# Set the default figure background color to white
plt.rcParams['figure.facecolor'] = 'white'

# Initialize JavaScript visualization in Jupyter notebooks
shap.initjs()

def get_basic_shap(model, X_train, X_test):
    
    # Create a SHAP explainer object
    explainer = shap.Explainer(model, X_train)

    # Calculate SHAP values for the test set
    shap_values = explainer(X_test)
    
    return shap_values

shap_values = get_basic_shap(logreg_noise['model'], logreg_noise['X_train'], logreg_noise['X_test'])

pred_df = logreg_noise['pred_df']

shap_values = get_basic_shap(logreg_noise1['model'], logreg_noise1['X_train'], logreg_noise1['X_test'])

pred_df = logreg_noise1['pred_df']

# Assuming 'overallStatus' is your target variable and the rest are features
X = ct_rawr.drop('overallStatus', axis=1)
y = ct_rawr['overallStatus']
xgb_noise = XGB_model(X,y)

X1 = ct_rawn.drop('overallStatus', axis=1)
y1 = ct_rawn['overallStatus']
xgb_noise1 = XGB_model(X1,y1)

indiv = xgb_noise
shap_values = get_basic_shap(indiv['model'], indiv['X_train'], indiv['X_test'])

pred_df = indiv['pred_df']

indiv = xgb_noise1
shap_values = get_basic_shap(indiv['model'], indiv['X_train'], indiv['X_test'])

pred_df = indiv['pred_df']

# Filter main_col to only include items that are in the columns of ct_rawr
org_col = [item for item in org_col if item in ct_rawr.columns]
lead_col = [item for item in lead_col if item in ct_rawr.columns]
continent_col = [item for item in continent_col if item in ct_rawr.columns]
country_col = [item for item in country_col if item in ct_rawr.columns]
mesh_col = [item for item in mesh_col if item in ct_rawr.columns]
bb_mesh = [item for item in bb_mesh if item in ct_rawr.columns]
phases_col = [item for item in phases_col if item in ct_rawr.columns]
enrollment_col_r = [item for item in enrollment_col_r if item in ct_rawr.columns]
main_col = [item for item in main_col if item in ct_rawr.columns]

import matplotlib.pyplot as plt
import numpy as np

def plot_feature_importance(feature_importance, top_n=30, threshold=0.5, var='main', model='Rare Diseases'):
    """
    Plots a horizontal bar plot of the mean_abs_shap_value_test column from the feature_importance DataFrame.
    
    Parameters:
    feature_importance (pd.DataFrame): DataFrame containing the feature importance values.
    top_n (int): Number of top features to display.
    threshold (float): Maximum value for the x-axis.
    var (str): Variable type.
    model (str): Model name.
    """
    # Sort the DataFrame by mean_abs_shap_value_test in descending order
    feature_importance = feature_importance.sort_values(by='mean_abs_shap_value_test', ascending=False)
    feature_importance = feature_importance.head(top_n)
    
    # Create the plot
    plt.figure(figsize=(8, 4))
    plt.barh(feature_importance.index, feature_importance['mean_abs_shap_value_test'])
    
    # Add labels and title
    plt.xlabel('Mean(|SHAP value|) average impact on model output magnitude')
    plt.ylabel('')
    plt.title('SHAP Feature Importance on ' + model + ' Test Set for ' + var + ' variables')
    
    # Set x-axis ticks
    if threshold <= 0.2:
        x_ticks = np.arange(0, threshold + 0.02, 0.02)
    else:
        x_ticks = np.arange(0, threshold + 0.05, 0.05)
    plt.xticks(x_ticks)
    
    # Add vertical grid lines
    plt.grid(axis='x', alpha=0.7)
    
    # Limit x range
    plt.xlim(0, threshold)
    
    # Invert y-axis to have the highest values at the top
    plt.gca().invert_yaxis()
    
    # Show the plot
    plt.show()

# Example usage:
# Assuming feature_importance is already defined
# plot_feature_importance(feature_importance, 10, 0.5)

shap_int = xgb_noise
shap_interpreter_r = ShapModelInterpreter(shap_int['model'])
feature_importance_r = shap_interpreter_r.fit_compute(shap_int['X_train'],shap_int['X_test'], shap_int['y_train'], shap_int['y_test'])
 
shap_int1 = xgb_noise1
shap_interpreter = ShapModelInterpreter(shap_int1['model'])
feature_importance = shap_interpreter.fit_compute(shap_int1['X_train'],shap_int1['X_test'], shap_int1['y_train'], shap_int1['y_test'])

feature_importance_r

plot_feature_importance(feature_importance_r,20,0.5)

ax2 = shap_interpreter_r.plot('summary', max_display=20)

ax1 = shap_interpreter.plot('summary',max_display=20)

plot_feature_importance(feature_importance,20,0.5,model='Common Diseases')

main_colx =['oversighthasdmc',
 'isfdaregulateddrug',
 'nb_primaryOutcomes',
 'nb_secondaryOutcomes',
 'original_enrollment_count',
 'conditions_count',
 'condition_mesh_count',
 'intervention_mesh_count',
 'primaryOutcomesSame',
 'secondaryOutcomesSame',
 'orgFullNameSame',
 'responsiblePartySame',
 'leadSponsorSame',
 'PRINCIPAL_INVESTIGATOR_responsibleParty',
 'SPONSOR_INVESTIGATOR_responsibleParty']
main_col_dict = {'oversighthasdmc':'Data Monitoring Committee',
 'isfdaregulateddrug':'FDA Regulated Drug',
 'nb_primaryOutcomes':'Number of Primary Outcomes',
 'nb_secondaryOutcomes':'Number of Secondary Outcomes',
 'original_enrollment_count':'Theorical Enrollment Count',
 'conditions_count':'Number of Conditions Keywords',
 'condition_mesh_count':'Number of MeSH Terms for Conditions',
 'intervention_mesh_count':'Number of MeSH Terms for Interventions',
 'primaryOutcomesSame':'No modification in Primary Outcomes',
 'secondaryOutcomesSame':'No modification in Secondary Outcomes',
 'orgFullNameSame':'No change of study organization',
 'responsiblePartySame':'No change of responsible party',
 'leadSponsorSame':'No change of lead sponsor',
 'PRINCIPAL_INVESTIGATOR_responsibleParty':'PI as responsible party',
 'SPONSOR_INVESTIGATOR_responsibleParty':'Sponsor-PI as responsible party'}

## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=main_col_dict,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(main_col_dict).values
shap_interpreter_r.column_names = [main_col_dict.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=main_col_dict,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(main_col_dict).values
shap_interpreter.column_names = [main_col_dict.get(col, col) for col in shap_interpreter.column_names]

## replace the main_colx with the dict in shap_interpreter_r.X_train
main_colx = [main_col_dict.get(item,item) for item in main_colx]
#rename columns names with the dict
ct_rawr = ct_rawr.rename(columns=main_col_dict)
ct_rawn = ct_rawn.rename(columns=main_col_dict)

plot_feature_importance(feature_importance_r.loc[main_colx],12,0.5,var='General Protocol')

ax1 = shap_interpreter_r.plot('summary', target_columns=main_colx,max_display=12)

plot_feature_importance(feature_importance.loc[main_colx],12,0.5,var='General Protocol',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=main_colx,max_display=12)

phases_col=['PHASE1','PHASE3']
xxx=exp_drug_col+phases_col

xxx_dict = {'experimental_only':'Experimental Drug Alone',
    'exp_against_placebo':'Experimental Drug against Placebo',
    'exp_against_active_comp':'Experimental Drug against Active Comparator',
    'PHASE1':'Phase 1',
    'PHASE3':'Phase 3'}

dictionnari = xxx_dict
## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter_r.column_names = [dictionnari.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter.column_names = [dictionnari.get(col, col) for col in shap_interpreter.column_names]

## replace the main_colx with the dict in shap_interpreter_r.X_train
xxx = [dictionnari.get(item,item) for item in dictionnari]
ct_rawr = ct_rawr.rename(columns=dictionnari)
ct_rawn = ct_rawn.rename(columns=dictionnari)

plot_feature_importance(feature_importance_r.loc[xxx],12,0.12,var='Multi-phases and experimental',model='Rare Diseases')

ax1 = shap_interpreter_r.plot('summary', target_columns=xxx,max_display=50)

plot_feature_importance(feature_importance.loc[xxx],12,0.12,var='Multi-phases and experimental',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=xxx,max_display=50)

elig_dict = {'interventionModel_CROSSOVER' : 'Cross-over Model',
 'interventionModel_PARALLEL':'Parallel Model',
 'interventionModel_SEQUENTIAL':'Sequential Model',
 'interventionModel_SINGLE_GROUP':'Single Group Model',
 'FEMALE_ONLY':'Female Only',
 'healthyVolunteers':'Healthy Volunteers eligible',
 'include_elderly':'Elderly eligible',
 'include_children':'Children eligible',
 'ci_pregnancy':'Pregnancy exclusion criteria',
 'arm_group_count':'count of arm groups','is_randomized':'Randomization'}
 
 


dictionnari = elig_dict
## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter_r.column_names = [dictionnari.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter.column_names = [dictionnari.get(col, col) for col in shap_interpreter.column_names]

elig_col = list(elig_dict.values())
elig_col
ct_rawr = ct_rawr.rename(columns=dictionnari)
ct_rawn = ct_rawn.rename(columns=dictionnari)

plot_feature_importance(feature_importance_r.loc[elig_col],10,0.22,var='Models and Eligibility',model='Rare Diseases')

ax1 = shap_interpreter_r.plot('summary', target_columns=elig_col,max_display=10)

plot_feature_importance(feature_importance.loc[elig_col],10,0.22,var='Models and Eligibility',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=elig_col,max_display=10)

continent_dict = {'continent_Australia and New Zealand': 'Australia and New Zealand',
 'continent_Central America': 'Central America',
 'continent_Central Asia': 'Central Asia',
 'continent_Eastern Africa': 'Eastern Africa',
 'continent_Eastern Asia': 'Eastern Asia',
 'continent_Eastern Europe': 'Eastern Europe',
 'continent_Northern Africa': 'Northern Africa',
 'continent_Northern America': 'Northern America',
 'continent_Northern Europe': 'Northern Europe',
 'continent_South America': 'South America',
 'continent_South-Eastern Asia': 'South-Eastern Asia',
 'continent_Southern Africa': 'Southern Africa',
 'continent_Southern Asia': 'Southern Asia',
 'continent_Southern Europe': 'Southern Europe',
 'continent_Western Africa': 'Western Africa',
 'continent_Western Asia':  'Western Asia',
 'total_centers'  : 'Count of Centers',
 'multi_city' : 'Centers in dinstinct Cities',
 'multi_country': 'Centers in distinct Countries',
 'multi_subregion':'Centers in distinct Subregions'}


 
 


dictionnari = continent_dict
## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter_r.column_names = [dictionnari.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter.column_names = [dictionnari.get(col, col) for col in shap_interpreter.column_names]

continent_col = list(dictionnari.values())
continent_col
ct_rawr = ct_rawr.rename(columns=dictionnari)
ct_rawn = ct_rawn.rename(columns=dictionnari)

plot_feature_importance(feature_importance_r.loc[continent_col],8,0.45,var='Location',model='Rare Diseases')

ax1 = shap_interpreter_r.plot('summary', target_columns=continent_col,max_display=8)

plot_feature_importance(feature_importance.loc[continent_col],8,0.45,var='Location',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=continent_col,max_display=8)

entity_dict = {'org_study_class_HOSPITAL':'Medical org study class',
 'org_study_class_INDUSTRY':'Industry org study class',
 'org_study_class_NETWORK':'Network org study class',
 'org_study_class_NIH':'NIH org study class',
 'org_study_class_UNIVERSITY':'Academic org study class',
 'sponsor_class_HOSPITAL' : 'Medical Lead Sponsor class',
 'sponsor_class_INDUSTRY':'Industry Lead Sponsor class',
 'sponsor_class_NETWORK':'Network Lead Sponsor class',
 'sponsor_class_NIH':'NIH Lead Sponsor class',
 'sponsor_class_UNIVERSITY':'Academic Lead Sponsor class',
 'number_of_collaborators':'Number of Collaborators',
 'industry_ratio':'Percentage of Industrial Collaborators',}


 
 


dictionnari = entity_dict
## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter_r.column_names = [dictionnari.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter.column_names = [dictionnari.get(col, col) for col in shap_interpreter.column_names]

entity_col = list(dictionnari.values())
entity_col
ct_rawr = ct_rawr.rename(columns=dictionnari)
ct_rawn = ct_rawn.rename(columns=dictionnari)

plot_feature_importance(feature_importance_r.loc[entity_col],8,0.18,var='Entity',model='Rare Diseases')

ax1 = shap_interpreter_r.plot('summary', target_columns=entity_col,max_display=8)

plot_feature_importance(feature_importance.loc[entity_col],8,0.18,var='Entity',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=entity_col,max_display=8)

bb_dict = {'bb_mesh:Blood and Lymph ': 'Blood and Lymph Conditions',
 'bb_mesh:Digestive System ':'Digestive System Diseases',
 'bb_mesh:Disease at or Before Birth':'Diseases and Abnormalities at or Before Birth',
 'bb_mesh:Gland and Hormone Related ':'Gland and Hormone Related Diseases',
 'bb_mesh:Heart and Blood ': 'Heart and Blood Diseases',
 'bb_mesh:Immune System ':'Immune System Diseases',
 'bb_mesh:Infections':'Infections',
 'bb_mesh:Musculoskeletal ':'Musculoskeletal Diseases',
 'bb_mesh:Neoplasms':'Neoplasms',
 'bb_mesh:Nervous System ':'Nervous System Diseases',
 'bb_mesh:Nutritional and Metabolic ':'Nutritional and Metabolic Diseases',
 'bb_mesh:Respiratory Tract':'Respiratory Tract (Lung and Bronchial) Diseases',
 'bb_mesh:Skin and Connective Tissue ':'Skin and Connective Tissue Diseases',
 'bb_mesh:Urinary Tract, Sexual Organs, and Pregnancy ':'Urinary Tract, Sexual Organs, and Pregnancy Conditions'}


 
 


dictionnari = bb_dict
## rename the indexes in feature_importance_r with the dict 
feature_importance_r.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter_r.X_test.columns = shap_interpreter_r.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter_r.column_names = [dictionnari.get(col, col) for col in shap_interpreter_r.column_names]

## rename the indexes in feature_importance_r with the dict 
feature_importance.rename(index=dictionnari,inplace=True)
## replace the columns names with the dict in shap_interpreter_r.X_test 
shap_interpreter.X_test.columns = shap_interpreter.X_test.columns.to_series().replace(dictionnari).values
shap_interpreter.column_names = [dictionnari.get(col, col) for col in shap_interpreter.column_names]

bb_mesh = list(dictionnari.values())
bb_mesh
ct_rawr = ct_rawr.rename(columns=dictionnari)
ct_rawn = ct_rawn.rename(columns=dictionnari)

plot_feature_importance(feature_importance_r.loc[bb_mesh],20,0.42,var='Main MeSH Terms',model='Rare Diseases')

ax1 = shap_interpreter_r.plot('summary', target_columns=bb_mesh,max_display=20)

plot_feature_importance(feature_importance.loc[bb_mesh],20,0.42,var='Main MeSH Terms',model='Common Diseases')

ax1 = shap_interpreter.plot('summary', target_columns=bb_mesh,max_display=20)

ax1 = shap_interpreter_r.plot('summary', target_columns=org_col,max_display=7)
ax2 = shap_interpreter_r.plot('importance', target_columns=org_col,max_display=7)

ax1 = shap_interpreter_r.plot('summary', target_columns=lead_col,max_display=7)
ax2 = shap_interpreter_r.plot('importance',target_columns=lead_col,max_display=7)

xgb_noise['name'] = 'XGB RD'
logreg_noise['name'] = 'LR RD'
xgb_noise1['name'] = 'XGB CD'
logreg_noise1['name'] = 'LR1 CD'

# Function to format metrics
def format_metric_with_std(row, metric):
    return f"{row[metric]:.3f} ± {row[metric + ' Std']:.2f}"

# Step 1: Create the DataFrame with Accuracy, Log-Loss, ROC-AUC, F1, PR-AUC, Recall, Precision, and their standard deviations

dictionaries = [logreg_noise['CV_scoring'], xgb_noise['CV_scoring'],logreg_noise1['CV_scoring'], xgb_noise1['CV_scoring']]
dictionary_name = ['logreg RD',  'xgb RD','logreg CD',  'xgb CD']
metrics = ['precision','f1','auc','pr_auc','accuracy', 'recall', 'log_loss']
cv_data = {
    'Model': dictionary_name  # ['logreg', 'xgb']
}
# Dynamically populate the DataFrame columns for each metric and its standard deviation
for metric in metrics:
    cv_data[metric] = [d[d['Metric'] == metric]['Mean Score'].values[0] for d in dictionaries]
    cv_data[f'{metric} Std'] = [d[d['Metric'] == metric]['Standard Deviation'].values[0] for d in dictionaries]

# Create the DataFrame from the constructed dictionary
df_models_cv = pd.DataFrame(cv_data)


# Step 2: Plot function to include all metrics with error bars
def plot_selected_scores_seaborn():
    # Define the metrics to plot
    metrics = ['precision','f1','auc','pr_auc','accuracy', 'recall', 'log_loss']
    
    # Create standard deviation metric names
    std_metrics = [f'{metric} Std' for metric in metrics]
    
    # Melt the DataFrame to make it suitable for Seaborn's barplot
    df_melted = df_models_cv.melt(id_vars="Model", value_vars=metrics, var_name="Score", value_name="Value")
    df_melted_std = df_models_cv.melt(id_vars="Model", value_vars=std_metrics, var_name="Score Std", value_name="Std")
    
    # Merge the melted DataFrames to have both Value and Std in the same DataFrame
    df_melted['Std'] = df_melted_std['Std']
    
    # Set the figure size
    plt.figure(figsize=(14, 8))
    
    # Create a bar plot with error bars
    barplot = sns.barplot(x='Model', y='Value', hue='Score', data=df_melted, ci=None)
    
    # Add error bars for each bar
    for i, bar in enumerate(barplot.patches):
        if bar.get_height() > 0:
            # Determine the corresponding model and metric (score)
            model = df_melted.iloc[i // len(metrics)]['Model']
            score = df_melted.iloc[i // len(metrics)]['Score']
            # Retrieve the corresponding standard deviation
            std = df_melted[(df_melted['Model'] == model) & (df_melted['Score'] == score)]['Std'].values[0]
            # Add the error bar
            barplot.errorbar(bar.get_x() + bar.get_width() / 2, bar.get_height(), yerr=std, fmt='none', c='black', capsize=5)

    # Annotate numerical values on top of each bar
    for p in barplot.patches:
        if p.get_height() > 0:
            barplot.annotate(format(p.get_height(), '.2f'), 
                            (p.get_x() + p.get_width() / 2., p.get_height() -.1), 
                            ha='center', va='bottom', 
                            xytext=(0, 10), 
                            textcoords='offset points')
    
    # Final touches for the plot
    plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
    plt.xlabel('Model')  # X-axis label
    plt.ylabel('Scores')  # Y-axis label
    plt.ylim(0, 1)  # Y-axis limit (assuming scores are normalized between 0 and 1)
    plt.title('PHASE 2 Rare Diseases RD and Common Diseases CD \nModel Performance Metrics with Standard Deviations')  # Title
    plt.tight_layout()  # Adjust layout to fit elements properly
    plt.legend(loc='upper center', shadow=True, ncol=7)
    plt.show()  # Display the plot

# Call the function to generate and display the plot
plot_selected_scores_seaborn()
# Create a new DataFrame with formatted metrics
formatted_data = {
    'Model': df_models_cv['Model']
}
for metric in metrics:
    formatted_data[metric] = df_models_cv.apply(lambda row: format_metric_with_std(row, metric), axis=1)
df_formatted = pd.DataFrame(formatted_data)
df_formatted